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Boolean networks, widely used to model gene regulation, exhibit a phase transition between 
regimes in which small perturbations either die out or grow exponentially. We show and numerically 
verify that this phase transition in the dynamics can be mapped onto a static percolation problem 
which predicts the long-time average Hamming distance between perturbed and unperturbed orbits. 



Boolean networks have been a prominent tool for mod- 
eling gene regulation since their introduction by Kauff- 
man in 1969 [H, 0). In a Boolean network, each node is 
assigned a state, or 1, which is synchronously updated 
at discrete time steps according to a pre-assigned update 
function which depends on the states of that node's in- 
puts on the previous time step. When used to model 
gene regulatory networks, each node represents a gene, 
and the state of the node indicates whether or not the 
gene is being expressed. Kauffman's original considered 
random networks and update functions in which each of 
the N nodes has K input links from randomly chosen 
nodes (the N-K model). Kauffman found numerically 
that when the in-degree K crosses a critical value, there 
is a transition between a stable phase, in which small per- 
turbations die out, to an unstable phase, in which small 
perturbations grow and become macroscopic. 

A derivation of the critical in-degree was given by 
Derrida and Pomcau for annealed N-K networks [3|. 
Here "annealed" means that the network edges and up- 
date functions are randomly redrawn between time steps. 
They hypothesized that for large networks the stability 
properties of the annealed system are similar to those of 
the original frozen (non-annealed) system. This hypoth- 
esis is well-supported by numerical experiments 0, Hj], 
and we refer to it as the "annealed approximation." Re- 
cent work [f| has extended this approach by using a par- 
tial randomization, in which only the update functions 
(but not the network topology) are randomly generated 
at each time step. In contrast with the annealed approx- 
imation, this "semi-annealed" approximation describes 
the dynamics on a fixed network which may have non- 
trivial topological features such as edge assortativity Q , 
motifs @ , and community structure jl] . The only neces- 
sary assumption is that the network is locally treelike (it 
cannot have many short loops) [i| . 

Some recent papers have derived stability properties 
of Boolean networks without annealing [HI Ej. These 
papers are complementary to ours in the following sense. 
Although rigorous, their results only apply to the ensem- 
ble average of random networks with restrictions on their 
network topology and/or update functions. In contrast, 
because our results rely on the semi-annealed approxima- 
tion, they can model the dynamics of a specific network. 

Here, using our semi-annealed approach, we map the 



dynamical problem of stability on a Boolean network onto 
the static problem of network percolation in the N — > oo 
limit. Previous authors have discussed the percolation 
properties of the "frozen component" of N-K networks 
13-14 1, and others have used percolation to discuss the 
stability of N-K lattices [TBI, [l6|. In contrast, we show 
that a dynamic quantity, the long-time average Ham- 
ming distance between two initially close trajectories on 
a Boolean network, can be mapped onto the size of the 
giant out-component in a percolation problem. We will 
illustrate this map in three different contexts. First, 
we consider the well-known annealed approximation and 
map it onto percolation in the configuration model 
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Second, we give a similar map from the semi-annealed 
approximation [f| to weighted site percolation [l8|. Fi- 
nally, we treat a more general class of update functions 
by mapping to a correlated bond percolation problem. 

Model: A Boolean network is a directed network of N 
nodes, in which each node i is assigned a state, Xi(t) = 
or Xi(t) = 1, at each discrete time step t. We denote the 
in- and out-degrees of node i by d™ and d° ut and the set 
of inputs to node i by Ji. A Boolean function or "truth 
table" Fi, fixed in time, updates the state of each node i 
at each time step t, Xi(t) = Fi {{xj(t — 1) : j € Ji})- 

In the literature, the truth tables Fi are usually gener- 
ated randomly (e.g., Q)- For each combination of input 
states to node i, the value of Fi is assigned to be 1 with 
probability p or with probability 1 — p, where p is the 
"bias probability." Below, as in |5j, we will consider the 
more general case where each Fi is generated with a dif- 
ferent bias pi assigned to each node i. Later, we will also 
consider the case of "canalizing" functions, in which one 
input acts as a master switch for the truth table. That 
is, input j to node i is canalizing if there is a state of 
Xj which completely determines the value of Fi indepen- 
dent of the other inputs to i. (When Xj is not equal to 
its canalizing value, F{ depends on the states of its other 
inputs.) Canalizing functions are thought to be common 
in real gene networks 3, 2fj| . 

Consider two trajectories, x(t) and x(t), which evolve 
on the same Boolean network. The initial conditions x(0) 
and x(0) differ only on a small randomly chosen fraction 
e of nodes. We say that a node i is "damaged" at time 
t if Xi(t) ^ Xi(t), and our goal is to predict the extent 
of the damage at long times. Let yi be the fraction of 
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time that node i is damaged on an orbit of length T as 
T — > oo. The normalized long-time average Hamming 
distance Y = (i/j), < Y < 1, is used as the order 
parameter for the stability phase transition. The average 
(•} is taken over all nodes i, then over all initial conditions 
which differ on a fraction e of the nodes. 

Analytic Results: First we treat the annealed approx- 
imation for random networks [2l|. We assume that the 
truth tables are randomly generated with a bias which 
depends only on degree. Let Pjk be the probability that 
a node has j inputs and k outputs, and let the bias of 
such a node be pjk- We define the sensitivity [Hf to be 
qjk = 2pjk (1 — Pjk) £ [0, 1], which can be interpreted as 
the probability that a node with j inputs and k outputs 
will become damaged at time t if at least one of its inputs 
is damaged at time t — 1. 

In the annealed approximation, Y can be predicted an- 
alytically using a method derived in Q and 0, which 
can be explained as follows. Let z denote the average 
degree of the network, i.e. z = J2 3 ,k3 p jk = Ej,fc k-P/fci 
and let E denote the probability that a randomly se- 
lected edge originates from a damaged node. A ran- 
domly selected edge originates from a node with j inputs 
and k outputs with probability p ' k , and such a node 
will become damaged with probability qjk if it has at 
least one damaged input, which occurs with probability 
1 - (1 - Ey. Therefore, 



degrees of nodes and edges, F (w) = ^ ■ k PjkqjkW^ and 
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Y = J2 P ^x [l-(l-£) J 



(1) 



In the stable regime, these equations only have the trivial 
solution E = and Y = 0, but there will be a nonzero 
solution in the unstable regime (23j . 

We now show that Eq. (TTJ) can be mapped onto the 
generating function formalism for treating weighted site 
percolation in directed configuration-model networks, as 
developed in [T^] and (24J. In this model, each node is 
deleted with some probability which depends only on its 
degree. The resulting ensemble of site-deleted networks 
exhibits a percolation phase transition, above which there 
is a macroscopic connected component or "giant compo- 
nent." This giant component contains a core of mutually 
path-connected nodes called the giant strongly connected 
component (GSCC); this, along with all the nodes which 
can be reached from it, is called the giant out-component 
(GOUT). In our map, we will identify the probability 
that a node is not deleted with the sensitivity, writing 
qjk for the probability that a node with j inputs and fc 
outputs is undeleted. With this identification, we will 
show that Y maps onto the expected fraction of nodes in 
GOUT, which we denote S. 

It is shown in [l?} and [24[ that S can be found as 
follows. First, define generating functions for the in- 



Fi(w) = J2j k~~r~ < li k ' w ' > ■ Next, let u be the probabil- 
ity that a randomly selected edge is not in GOUT. The 
authors show through diagrammatic expansion that 



u = l-F l {l)+F 1 {u), 
S = F (l) - F (u). 



(2) 



We note that the substitutions E = 1 — u and Y = S map 
Eq. |T]) onto Eq. Therefore, the phase transition be- 
tween dynamical stability and instability in this ensemble 
of random Boolean networks is equivalent to the static 
percolation phase transition on the same ensemble. 

Our second result is a more general derivation of the 
same correspondence, using the framework of This 
framework applies to a specific locally treelike network in 
which each node i can have its own arbitrarily chosen bias 
Pi, with an associated sensitivity qi = 2pi(l —pi). Again, 
we will identify the sensitivity qi with a site nondeletion 
probability and map the Hamming distance, Y, onto the 
size of GOUT, S. We begin by writing an analogue of Eq. 
((T|) for a specific node in a semi-annealed, locally treelike 
Boolean network, 



Vi = 1% 
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(3) 



This is the long-time limit of a damage-spreading equa- 
tion derived in [5j, which noted that i will become dam- 
aged with probability qi if at least one of its inputs is 
damaged. The assumption that the network is locally 
treelike is necessary because all the probabilities in the 
product are treated as independent. 

Reference [3| derives a similar condition for site per- 
colation on locally treelike directed networks in which the 
probability that each node is not deleted is qi. It defines 
rji as the fraction of site-deleted networks for which node 
i is not in GOUT, and it shows that 
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because a node is not in GOUT when it is either deleted 
or has no inputs from GOUT. We note that substituting 
yi = 1 — T]i maps Eq. (|3|) onto Eq. ((4|. Because Y = (iji) 
and S = (1 — rji), this map also yields Y = S. For 5, 
the average (•) is first taken over all nodes i, then over 
all node deletion trials. 

We now introduce a third case, in which we consider 
Boolean networks with canalizing functions. The method 
used for our previous results can be extended to canal- 
izing functions, but because the truth table elements in 
a canalizing function are not generated independently, 
we need to consider a new type of percolation problem 
which we call correlated bond percolation. Instead of 
typical bond percolation, in which each bond is occupied 
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FIG. 1 (color online). The ensemble averages of Y, S, and 
T (taken over 20 networks) versus the average degree z, for 
three families of networks. The three families of networks are 
assortative (left), neutral (middle), and disassortative (right). 



or deleted independently, we consider joint probabilities 
where the deletion of two bonds may be correlated if they 
are inputs to the same node. 

Here we describe a correlated bond percolation prob- 
lem that corresponds to a Boolean network whose truth 
tables each have one canalizing input but are otherwise 
generated randomly. That is, for each node i, there is a 
canalizing input Cj, and all the rows of the truth table 
on which x Ci assumes its canalizing value have the same 
constant output; but the outputs of the other rows are 
randomly generated with a probability bias pi. To be- 
gin, we imagine that the system is equally likely to be 
in any of its states. As we will show, it is then formally 
possible to obtain equations describing damage spread- 
ing in closed form. Based on our numerical results, we 
conjecture that these equations can be used to predict 
damage spreading in a large class of Boolean networks 
with frozen truth tables. 

Working under the supposition that all system states 
are equally probable, we now derive an expression for y^. 
Let ri denote the "activity" of Ci on i 22j , defined as the 
fraction of states in which i will become damaged if Cj 
becomes damaged. If Ci is not damaged, it may be in 
cither the canalizing or non-canalizing state, each with 
probability i. In the first case it is impossible for i to 
become damaged, while the second case is equivalent to 
Eq. ©. Therefore, 



FIG. 2 (color online). Y versus S for individual neutrally 
assortative networks. 



Vi = ny Ci + -qi (1 -y Ci ) 
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where J[ = Ji — \ci\ and qi is the sensitivity of the half 
of the truth table where x Ci is not in its canalizing state. 



It can be shown that this is equivalent to 



vi = i-n 



Vci + -jq% ]J Vj> 



(0) 



where r\i = X—yi. This corresponds to a correlated bond 
percolation problem in which one of the following three 
things may occur. With probability 1 — ri, all edges to 
i are deleted; with probability ri — \qi, all of i's edges 
are deleted except for the edge from cf, and otherwise no 
input edges arc deleted. Note that it is straightforward 
to describe the case where only some of the nodes have a 
canalizing input by using Eqs. (JSJIHJ) for those nodes and 
Eqs. (|3"1I4"|) for the others. 

Numerical Results: We begin with the map described 
by Eqs. (J3]l4j> , since it is more general than Eqs. (JT][2|). 
We compare the long-time average Hamming distance Y 
to the size of the giant out-component S for particular 
networks. We also compare both Y and S to the theo- 
retical prediction given by the solution to Eq. ([3]), which 
we denote T. 

Our algorithm is as follows. First we create a 
configuration-model network with N = 10 5 nodes. The 
data in the figures were obtained using networks with 
Poisson-distributcd in-degrees and scale- free out-degrees; 
we have also tested other degree distributions and found 
similar results. If desired, we then enhance interesting 
topological features such as assortativity or feedforward 
loops using the same algorithms as in Next, we assign 
each node a bias Pi. These may be distributed randomly, 
or, if we wish to encourage (impede) instability on the 
network, we distribute them so that the nodal average 
(qia™d° nt ) is maximized (minimized) For the data 
in the figures, the biases pi were distributed randomly so 
that the sensitivities qi form a uniform distribution on 
the interval [.3, .5]. We choose random initial conditions 
for x, and a randomly selected fraction e = .01 of the 
nodes arc flipped for the initial conditions of x. 
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FIG. 3 (color online), (a) Linear and (b) log-log scatterplots 
of Y versus S for data generated in the same way as that of 
Fig. 2, except that now we average over the quenched disorder 
in the truth tables as described in the text, (c) Linear and (d) 
log-log scatterplots of Y and S versus T for the same data, 
sampling alternate points for visibility. 



To find Y, we time-evolve the system and average 
\xi(t) — Xi(t)\ between t = 900 and t = 1000, averaging 
over 100 initial conditions. The theoretical prediction is 
found by iterating Eq. §S§ until it converges to a solution 
y, then taking T = (yj). Finding S is less straightfor- 
ward, because a typical percolation problem is only guar- 
anteed to have a single, well-defined giant out-component 
in the N — > oo limit. For reasons discussed in the online 
Supplemental Material, we choose the following proce- 
dure. We delete each node i with probability 1 — % and 
find any strongly connected components (SCCs) in the 
resulting network, where we define an SCC to be a mu- 
tually path-connected set of nodes containing at least one 
loop. We define S to be the fraction of nodes which can 
be reached from at least one SCC, averaged over the en- 
semble of deletion trials. We average 10 3 deletion trials 
per network. We find that the numerical uncertainty in 
our measured values of T, Y , and S for each point in Figs. 
1-4 is smaller than the point size; see the Supplemental 
Material for details. 

Figure 1 illustrates the relationship between Y, S, and 
T for networks generated in this way. We see that Y 
and S have the same average values on the ensemble of 
random networks with given average degree z. However, 
in Fig. 2, we see that the prediction Y — S sometimes 
fails for individual networks, especially near the phase 
transition. The deviations in Fig. 2 are primarily caused 
by the quenched disorder in the truth tables, which may 
cause orbits to fall onto attractors which visit only a small 
fraction of the state space (and so may deviate from the 
semi-annealed approximation) . 

In Fig. 3, we have averaged over this quenched disor- 
der by choosing truth tables from an ensemble of closely 
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FIG. 4 (color online). Scatterplots of Y versus S for networks 
in which each node has one canalizing input, using Eqs. (|5I6|I . 

related frozen truth tables (but not networks) as follows. 
Before wc time-evolve each new pair of initial conditions, 
we perform a set of exchanges on the truth tables. For 
each edge j —> i, with probability |, we exchange Xj = 
and Xj = 1 on the truth table for i. We note that 
there are two major differences between this and the 
semi-annealed approximation. In the latter, the truth ta- 
bles are changed during the dynamics, whereas here they 
are only changed before each new dynamical trial. Sec- 
ond, whereas the semi-annealed approximation treats all 
inputs interchangeably, this procedure preserves input- 
specific information (such as whether an input is canaliz- 
ing). In Fig. 3, we see that this procedure yields excellent 
agreement between Y , S, and T for individual networks 
well above the transition. Near the transition and below 
it, finite-size effects still cause S (and, to a lesser extent, 
Y) to deviate slightly from the prediction T. These ef- 
fects are discussed in the Supplemental Material. 

In Fig. 4, we perform the same numerical experiment 
for the case in which each node has one canalizing input. 
We find that Y, S, and T agree for individual networks 
when we use the map between Eqs. ([5]) and (|6]), but the 
map between Eqs. ([3]) and ((H) fails for this case, indicat- 
ing that we retain significant input-specific information 
about the dynamics when we average over the quenched 
disorder in the truth tables. 

Discussion: We have presented evidence that the sta- 
bility of a Boolean network can be understood in terms of 
a related percolation problem on that network. This rela- 
tionship may be helpful in understanding the stability of 
systems modeled by Boolean networks, such as gene reg- 
ulatory networks and neural networks. Two previously- 
studied cases (the annealed and semi-annealed approxi- 
mations) map onto known results for percolation, and a 
case of biological interest (canalizing truth tables) maps 
onto a novel percolation problem. These maps are valid 
for the typical cases in the literature (large, locally tree- 
like networks with random or canalizing truth tables), 
but have the advantage of applying to specific networks 
rather than ensembles of random networks. Numerical 
experiments show excellent agreement with our analysis 
when averaged over a family of quenched truth tables. 
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